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ABSTRACT 

This  paper  discusses  the  efficient  numerical  treatment  of  optimal  control  problems  gov¬ 
erned  by  elliptic  PDE’s  and  systems  of  elliptic  PDE’s,  where  the  control  is  finite  dimensional. 
Distributed  control  as  well  as  boundary  control  cases  axe  discussed.  The  main  characteristic 
of  the  new  methods  is  that  they  are  designed  to  solve  the  full  optimization  problem  directly, 
rather  than  accelerating  a  descent  method  by  an  efficient  multigrid  solver  for  the  equations 
involved.  The  methods  use  the  adjoint  state  in  order  to  achieve  efficient  smoother  and  a 
robust  coarsening  strategy.  The  main  idea  is  the  treatment  of  the  control  variables  on  ap¬ 
propriate  scales,  i.e.,  control  variables  that  correspond  to  smooth  functions  are  solved  for  on 
coarse  grids  depending  on  the  smoothness  of  these  functions.  Solution  of  the  control  prob¬ 
lems  is  achieved  with  the  cost  of  solving  the  constraint  equations  about  two  to  three  times 
(by  a  multigrid  solver).  Numerical  examples  demonstrate  the  effectiveness  of  the  method 
proposed  in  distributed  control  case,  pointwise  control  and  boundary  control  problems. 
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1  Introduction 


Computations  in  optimal  control  for  distributed  parameter  systems  demand  huge  com¬ 
putational  resources.  These  are  optimization  problems  with  constraints  which  are  par¬ 
tial  differential  equations.  This  paper  focuses  on  the  computational  aspects  of  optimal 
control  problems  governed  by  elliptic  partial  differential  equations  or  elliptic  systems, 
where  the  control  variables  are  of  finite  dimension. 

The  classical  optimal  control  algorithms  use  the  costate  (  a  Lagrange  multiplier)  in 
order  to  define  an  iteration  that  converges  to  the  minimum.  These  are  descent  type 
algorithms  which  involve  the  solution  of  the  state  and  costate  equations  per  iteration  of 
updating  the  control  variables.  The  state  and  costate  equations,  being  discretizations 
of  partial  differential  equations  are  quite  expensive  to  solve.  The  over  all  solution  of 
the  optimization  problem  becomes  therefore  very  expensive. 

The  need  for  more  efficient  methods  is  obvious.  One  way  of  improving  the  situation 
is  by  using  fast  solvers  ,  i.e.,  multigrid  methods,  for  the  state  and  costate  equations  at 
each  step  of  a  descent  algorithm.  This  may  result  in  a  significant  saving  in  compute 
tional  time,  but  will  not  affect  the  number  of  outer  iterations  for  updating  the  control 
unknowns’,  it  will  just  make  the  time  required  for  each  iteration  smaller. 

A  much  more  effective  solution  process  can  be  obtained  by  having  an  algorithm 
that  is  aimed  at  the  full  optimization  problem  directly,  rather  than  accelerating  a 
linear  equation  solver.  This  paper  is  devoted  to  this  approach.  Multigrid  algorithms 
that  were  designed  to  tackle  the  full  problem  rather  than  being  a  mere  fast  linear  solver 
have  been  developed  for  a  class  of  problems.  The  first  one  to  be  developed  was  the 
FAS  algorithm  for  nonlinear  equations  [Bl].  Recently  such  an  approach  was  successful 
in  treating  some  bifurcation  problems  [T2]  and  stability  calculation  [T2].  In  all  these 
cases  maximal  efficiency  has  been  obtained  as  a  result  of  treating  the  full  problem 
directly. 

The  high  efficiency  in  the  problems  treated  here  is  achieved  by  working  on  all 
unknowns  of  the  problems  (i.e.,  state,  costate  and  control  unknowns)  simultaneously 
where  scales  of  the  different  unknowns  are  taken  into  account.  The  main  rule  is  to 
treat  different  control  variables  on  different  grids  depending  on  their  scale  of  influence. 
This  means  that  a  variable  that  has  a  non  local  effect  on  the  solution  should  be  treated 
mainly  on  coarse  levels.  On  the  other  hand,  a  variable  that  has  a  non-smooth  effect 
on  the  solution  in  some  neighborhood,  will  be  relaxed  in  that  neighborhood  on  fine 
levels  as  well.  More  generally,  local  work  on  fine  grids  will  be  done  to  update  a  certain 
control  variable  only  at  the  vicinity  where  it  has  a  non  smooth  effect  on  the  solution 


while  work  on  coarser  levels  for  that  variable  will  be  in  larger  and  larger  neighborhoods 
of  that  region. 

As  is  done  in  the  classical  optimal  control  algorithms  a  costate  variable  is  introduced 
in  our  solution  process.  This  allows  the  construction  of  good  relaxation  schemes  in 
general,  and  is  even  more  important  when  coarsening  is  considered  (as  explained  in 
section  4).  In  the  relaxation  part  of  the  algorithm  we  distinguish  two  main  steps.  One 
is  designed  to  smooth  the  errors  in  the  state  and  costate  for  a  given  choice  of  the  control 
variables.  This  is  done  for  the  state  variable  using  the  Gauss-Seidel  (GS)  relaxations 
for  scalar  elliptic  equations  of  state  and  Distributed-Gauss-Seidel  (DGS)  for  systems  of 
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elliptic  equations  [B2].  The  equations  obtained  for  the  costate  are  very  similar  to  the 
state  equations  and  are  relaxed  using  essentially  the  same  relaxation  method  as  for  the 
state  equations.  The  other  step  of  the  relaxation  is  designed  for  updating  the  control 
variables.  This  step  of  the  relaxation  is  done  taking  into  account  the  smoothness  of 
the  change  introduced  in  the  state  and  costate  as  a  result  of  a  small  change  in  one  of 
the  control  variables.  That  is,  some  work  is  to  be  done  on  coarse  grids  only  and  some 
on  fine  levels  in  special  subdomains.  The  update  of  the  control  variables  is  done  in  the 
following  way.  A  small  change  is  introduced  in  them,  then  its  effect  on  the  state  and 
the  costate  is  computed  approximately.  The  amplitude  of  that  change  in  the  control 
is  then  chosen  so  as  to  reduce  the  functional  or  to  minimize  the  residuals  of  certain 
equations. 

The  algorithms  constructed  that  way  are  aimed  at  the  full  problem  directly  and 
need  no  iterations  within  iterations  as  the  standard  algorithms  require. 

Numerical  examples  show  that  solutions  to  the  levels  of  discretization  errors  are 
obtained  in  just  a  few  work  units,  where  a  work  unit  is  the  work  involved  in  relaxing 
the  state  equations.  The  complexity  of  solving  a  control  problem  is  about  twice  that  of 
solving  the  corresponding  constraint  equations,  which  is  essentially  optimal.  Examples 
for  distributed  control,  boundary  control  and  pointwise  control  are  given  and  all  show 
the  same  efficiency. 

2  The  Problem 

Let  V,W,W  and  H  be  Hilbert  spaces.  Let  A,B  and  C  be  operators  as  follows  A  :  V  -* 
W,  B  :U  -»  W,  C  :V  -»  K  and  let  Uad  C  U.  Consider  the  problem 

\\Cx  -  d||’  .  . 

Ai  =  Bu  +  f  K  } 

Here  A  is  assumed  to  be  an  elliptic  partial  differential  equation  (PDE)  or  an  elliptic 
system  of  PDE’s  with  smooth  coefficients  (see  [ADN]  ).  B  is  assumed  to  be  a  finite 
dimensional  operator,  that  is,  there  exist  functions  if/j  G  W  j  =  1 . .  such  that 

Bu=^2  (2-2) 

i=i 

We  assume  that  the  functions  ipj  are  smooth  except  for  a  finite  number  of  singularities. 
U  is  then  identified  with  IK.  The  adjoint  of  B,  B*  :  W  — *  IK  is  given  by 

B*<t>  -  (<  ^1,0  >w»---»<  'I’vt  >w)T  (2-3) 

The  following  equations  hold  at  the  minimum, 

Ax  =  Bu  +  / 

A*p  +  C*Cx  =  C'd  (2.4) 

B*p  =  0 
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where  B  and  B*  are  given  by  (2.2)  and  (2.3)  respectively.  Let  xq  be  a  solution  of 
Axo  =  f  and  d  =  d  —  Cxq.  Problem  (2.1  )  can  be  written  in  terms  of  x  =  £  -  *o  and 

it  as 


minuew.d  IIC1*  -  d\\2n 
Ax  =  Bu 

So  without  loss  of  generality  we  can  assume  that  /  =  0. 


(2.5) 


2.1  Descent  Methods 

Iterative  methods  for  solving  the  control  problem  are  described.  A  triple  (x,u,p)  will 
be  called  compatible  if  it  satisfies  the  equations 


Ax  -  Bu  =  0 
A0p  +  C*Cx  =  C*d 


(2.6) 


Given  an  approximate  solution  (s,  u)  of  (2.5),  an  improved  one  is  obtained  as  follows. 
Denote  by  E(x,u)  the  functional  DC's  -  d||2 .  For  a  compatible  triple  (x,u,p) 

E(x,  u)  =  -  <  Cx,Cx  >H  -2  <  u,B*p  >u  +  <  d,d  >n  (2.7) 

We  look  for  a  change  711  in  u  and  corresponding  changes  in  x,p  that  result  in  a  new 
compatible  triple,  and  for  which  the  value  of  E  is  smaller.  Let  (s  +  7s,  u  +  7 u,p  +  ip) 
be  a  new  compatible  triple.  Then 


E(  x  +  7  i,u  +  7«)  =  E(x,  tt)-f 

-27  <  Cx,Cx  >n  -72  <  Cx,C£  >n  (2.8) 

-27  <  B*p,u  >u  -27  <  B*p,u  >u  -272  <  B*p,  u>n 


Choosing  u  =  B*p,  and  £,p  satisfying  Ax  -  Bu  =  0  and  A*p  +  C*Cx  =  0  we  get 


E{x  +  7i,  u  +  7ft)  =  E(x,u) 

-27(<  Cx,Cx  >H  +  <  B*p,B*p  >u  +  <  B*p, u  >u ) 

-72(2  <  B*p,B*p  >u  +  <  Cx,Cx  >„) 

Among  all  possible  7,  the  choice 

_  1  <  Cx,Cx  >H  +  <  B*p,B*p  >u  +  <  B*p}u  >u 

7_  2  |jC5:||2  +2<B*p,B*p>u 

minimizes  the  functional  in  the  direction  chosen.  ThiB  can  be  a  basis  for  an  algorithm 
for  solving  the  optimal  control  problem. 

A  Descent  Algorithm 


(2-») 


(2.10) 
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Set  «o  =  O.po  =  0,7  =  1,  n  =  1 . 

While  7  >  e  Do 
Begin 

1.  Solve  for  x„  the  equation  Axn  -  Bun  =  0 

2.  Solve  for  pn  the  equation  A*pn  + C*(Cxn  — d)  =  0 

3.  Set  u  =  B*pn .  Solve  Ax  -  Bu  =  0,  A*p  +  C*Cx  =  0. 
(approximately)  for  ®,p  respectively. 

4.  Set  un+1  =u„  +  7u,  in+i=®n  +  7i,  pn+i  =  Pn  +  7P.  n  =  n+l 
with  7  given  by  (2 . 10) . 

End 


Approximate  Descent  Algorithm 

Set  uo  =  O.po  =  0,7  =  l,n  =  1. 

While  7  >  c  Do 
Begin 

1.  Starting  with  xn  =  ®n_x  relax  Axn  =  Bun,  yielding 
a  solution  satisfying  Axn  ~  Bun  =  r” 

2.  Starting  withp„=p„_i  relax  A*pn  +  C*Cxn  =  C*d, 
yielding  a  solution  satisfying  A*Pn  +  C*(Cxn  -  d)  = 

3.  Set  u  =  B*pn.  Solve  Ax  -  Bu  =  0,  A*p  +  C*Cx  =  0 
(approximately)  for  £,p  respectively. 

4.  Set  txn+i  =  ti„  +  7u,  ®Tl+i=®n  +  7**  Pn+i=Pn  +  n  =  n+l 
with  7  given  by  (2 . 10) . 

End 

e  is  a  small  enough  number,  depending  on  the  accuracy  required  for  the  solution. 
Note  that  if  (|r£(|  -»  0  and  ||r”||  — ►  0  then  the  above  algorithm  converges  to  the 
minimum  of  (2.5). 

The  expensive  steps  in  this  algorithm  are  the  ones  for  obtaining  z  and  p  by  ap¬ 
proximately  inverting  A  and  A *,  which  are  elliptic  operators  (or  systems).  Moreover, 
the  solution  of  the  elliptic  equations  governing  the  state  and  costate  needs  to  be  done 
many  times. 

Remark  Considering  the  form  of  the  necessary  conditions  one  may  be  tempted 
to  choose  7  based  on  minimizing  the  L%  residuals  of  the  third  equations  of  (2.4).  This 
leads  to 


<  B*p,B*p  > 

<  B*p,B*p  > 


(2.11) 
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which  may  not  be  a  good  choice  as  <  B*p,  B*p  >  may  vanish,  causing  the  iteration  to 
get  stuck  away  from  the  minimum. 

3  Naive  Multigrid  Approaches 

The  approximate-descent  algorithm  presented  in  the  previous  section  suggests  an  im¬ 
mediate  acceleration  procedure,  namely,  using  some  fast  solver  in  steps  1,2,  and  3. 
Since  the  form  of  the  operator  A*  is  very  similar  to  that  of  A  the  same  fast  solver  can 
be  used  in  all  the  steps.  The  operator  A  is  assumed  to  be  elliptic  and  therefore  calls 
for  multigrid  methods. 

Assume  that  we  are  given  V*,W*,W*  and  Hk  Hilbert  spaces  approximation  the 
original  V,  W,U  and  H.  Let  A*,  5*  and  Ck  be  operators  as  follows  Ak  :  V*  — ♦  W*, 
Bk'-Uk  -*  Wfc,  Ck  "Vk  — ►  'Hk-  Assume  also  that  interpolation  and  restriction  operators 

***  given;  Rkx  :  Vfc+i  -  Vtl  B*  :  Wfc+1  -  W*,  Rk  :  Uk+i  -» Uk , 
Pk  :  Vk  -*  V*+i,  Pk  :  Wfc  — »  Wfc+i,  Pk  :  Uk  -»  Uk+i-  The  superscript  k  in  these 
expressions  represent  a  level  of  discretization  where  k  =  1  correspond  to  the  coarsest 
level. 

Denote  by  MG(Ak,b,x,y,v\,V2)  a  multigrid  V(i/i,*/a)-cycle  for  solving  AkX  =  b  on 
level  k,  starting  with  initial  approximation  x  ending  with  y.  Assume  that  a  solution  of 
the  minimization  problem  is  needed  for  level  k  =  m.  A  possible  fast  algorithm  for  it, 
starting  with  initial  approximation  (zm,  um,pm)  is  denoted  by 


(xm,um,pm)  *-  MinMG(Am,Bm,Cm,xm,um,pm,<r) 


(3.12) 


and  given  by 


Set  u™  =  O.Pq1  =  0,  7  =l,n  =  l. 

While  7  >  e  Do 
Begin 

1.  Perform  a  cycle  MG(Am,flmu™,x™,z™+1, 17,1/2) 

2.  Perform  a  cycle  MG{A^,C^<T  -  C^Cx^^^+i, vuv2) 

3.  Perform  a  cycle  MG(4,5mOm,*”  1,*'*) 

4.  Perform  a  cycle  MG(A^,-C^Cm*”+1,p^,p^+1,»/i,i/2) 

6.  Set  <+i  =tC  +  7*C+i.  *”+i  =*«  +7*Sn.  Pn+i=Pn  +7#Cn. 
n  =  n  -|- 1  and  7  given  by  (2.10)  . 

End 

An  obvious  improvement  of  this  algorithm  is  to  start  with  a  good  initial  approxi¬ 
mation  for  the  optimization  problem  on  level  m.  This  can  be  achieved  easily  by  solving 
the  problem  first  on  level  m  -  1.  Applying  this  idea  recursively  one  arrives  at  the 
following  algorithm. 
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Set  =  0,®J  =  0,pj  =  0,k  =  1. 

While  k  <  m  Do 
Begin 

1.  Perform  a  cycle 

(®\u *,p*)  <-  MinMG(Ak,Bk,Ck,xk,uk,pk,dk ) 

2.  If  k  =  m  stop 

3.  ®*+1  =  n*®*,p*+1  =  n*p*,«*+1  =  p*uk,k  =  k  + 1. 

End 

where  II*,  11*  are  operators  analogous  to  P*  and  P*,  respectively,  but  of  possibly  a 
higher  order  (as  is  usual  in  multigrid  algorithms  for  elliptic  problems). 

Our  aim  is  to  develop  an  algorithm  in  which  the  total  computational  cost  of  solving 
the  control  problem  will  be  only  few  times  (2-3)  more  than  that  of  solving  for  the  state 
variable  x  given  that  u  is  known.  This  is  done  in  the  next  sections  where  a  multigrid 
approach  which  is  aimed  at  the  full  problem  directly  is  constructed. 

4  “One  Shot”  Multigrid  Methods 

In  this  section  we  describe  a  multigrid  method  that  is  aimed  at  the  full  optimization 
problem,  rather  than  serving  as  a  fast  solver  in  a  step  of  a  basic  optimization  algorithm, 
as  was  done  in  the  previous  section. 

The  construction  of  this  method  has  been  achieved  basically  by  following  the  next 
five  steps. 

1.  Distinction  between  smooth  versus  non-smooth  components  in  all  types  of  vari¬ 
ables  involved 

2.  Construction  of  a  basic  non-expensive  relaxation  for  the  full  problem. 

3.  Classification  of  rate  of  convergence  of  the  basic  relaxation  for  the  different  vari¬ 
ables  in  the  problem,  with  regard  to  their  smoothness. 

4.  Based  on  3,  determination  of  the  role  of  coarse  grids  in  accelerating  the  fine  grid 
convergence. 

5.  Construction  of  a  coarse  grid  problem  that  approximate  the  Pull  fine  grid  residual 
problem. 

Steps  1,2  and  3  are  discussed  in  detail  in  section  4.2,  the  rest  in  section  4.1.  The 
resulting  algorithm  involves  a  sequence  of  discrete  optimization  problems,  starting 
from  a  very  coarse  discretization  ending  with  a  fine  grid  one  on  which  a  solution  of  the 
problem  is  required.  The  optimization  problem  on  each  level  is  served  to  accelerate  the 
convergence  of  the  next  finer  level  problem.  The  iterative  process  for  solving  the  finest 
grid  problem  involves  a  relaxation  process  on  each  level  whose  effect  is  to  damp  out  the 
non-smooth  part  of  the  error,  i.e.,  the  part  of  the  error  which  cannot  be  represented 
on  the  next  coarse  grid  in  the  sequence.  This  together  with  a  proper  transferring  of 
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information  (residuals  and  current  line  grid  solution)  from  fine  to  coarse  levels  and 
back  to  fine  levels  results  in  a  very  high  efficiency.  A  full  description  of  the  method 
includes  a  relaxation  scheme  and  a  coarsening  strategy. 

4.1  Coarsening 

A  coarsening  scheme  requires  the  definition  of  coarse  grid  spaces  analogous  to  V,  W,  U,  H 
which  we  denote  by  the  same  letters  with  a  superscript  c.  Also  inter  grid  transfer 
P*,Pp,P»  and  Rx,Rpt  and  coarse  grid  operators  Ac,  Bc,  Cc  are  required.  Having  all 
this  we  can  define  the  coarse  grid  minimization  problem.  An  attempt  to  use 

min^ea^,  \\Cf{x*  +  Pxxc )  -  df  ||2 
AcXc  -  Bcuc  =  Rx(gf  -  Afx*  +  Bfu*)  '  ' 

as  an  approximation  for  the  error  on  the  fine  grid  fails.  The  reason  for  this  is  that 
this  coarse  grid  minimization  problem  does  not  preserve  the  fine  grid  minimum  once 
it  was  reached.  This  follows  immediately  from  the  necessary  conditions  of  the  coarse 
grid  problem  (4.1),  namely, 


AcXc  -  Bcuc  —  Rx(gl  ~  Afx*  +  Bfuf) 

A*epe  +  C*Ccxc  =  RvC){df  -  Cfx* )  (4.2) 

B*pc  =  0. 

The  coarse  grid  minimization  problem  formulated  above  is  a  problem  for  the  correc¬ 
tion  in  the  fine  grid  approximation.  Therefore,  upon  reaching  the  fine  grid  minimum 
the  coarse  grid  correction  xc,uc  should  be  zero.  This  however  is  not  the  case  since  the 
right  hand  side  of  the  second  equation  is  non-zero  in  general  at  the  minimum  of  the 
fine  grid  problem. 

In  order  to  define  a  correct  coarse  grid  problem  we  rewrite  the  fine  grid  minimization 
problem  as 


<  Cfxf’xf  >  -2  <  xf>4  >  "2  <  uf’9l  >  .  3x 

Afxt  -  Bfuf  =  gf  (4>3) 

where  g£  =  g,  g*  =  C*d  and  g*  =  0  on  the  fine  grid.  Note  that  the  problem  is  defined 
once  gltgl,9p  are  given. 

1.  Correction  Scheme.  The  coarse  grid  problem  for  the  correction  is  of  the  same 
form,  that  is, 


<  £e*e,*e  >  -2  <  xc,g‘  >  -2  <  uc,gcu  > 
AeXc  -  Beuc  =  g% 


(4.4) 
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where 


ffx  =  &x(g£  -  A/xc  +  Bfu f) 

9Cp  =  &r(gfP  ~  Alp*  -  C*cCcxf)  (4.5) 

gZ  =  Ru{gt  -  BfP*). 

The  necessary  conditions  for  this  coarse  grid  problem  are 

Acxc  -  Bcuc  =  g% 

A*PC  +  C*cCcxc  =  g%  (4.6) 

Kvc  =  aZ 

which  are  the  Correction  Scheme  (CS)  (see  [B2])  equations  for  the  fine  grid  necessary 
conditions.  In  particular  they  have  the  property  of  introducing  no  change  to  the  fine 
grid  problem  once  the  minimum  has  been  reached. 

Once  an  approximation  to  the  coarse  grid  problem  is  found  it  is  used  to  correct  the 
fine  grid  approximation  by 


x*  <—  X*  +  Pxxe 

uf  i-  u*  +  Puuc  (4.7) 

pf  *-pf  +  Pppc 

Note  that  this  coarse  grid  problem  depends  on  the  costate  p  obtained  after  the  fine 
grid  relaxation  process.  Without  using  it  the  coarsening  becomes  very  difficult  as  was 
explained  before. 

2.  Full  Approximation  Scheme.  The  coarse  grid  problem  for  quantities  that 
approximate  the  full  fine  grid  solution  rather  than  the  correction  for  it  have  similar 
form.  The  minimization  problem  is 


minu*ewcd  <  Cexc,xe  >  -2  <  x c,g*  >  -2  <  uc,gl  > 

Ac  xc  -  Bcuc  =  g% 

where 

g%  =  Mg£  -  Afxc  +  b}u*)  +  acrxx /  -  bc r^u* 
gcp  =  Rp(gl  -  A)pf  -  C*Cexf)  +  AIRppf  +  ctcj^xf 
SZ  =  £«($«  -  BfP*)  +  BlRpP*. 

The  operators  #*,  Rp,  R^  are  fine  to  coarse  grid  transfers  possibly  different  from  iZ»,  Rp,  Ru 
(see  [B2]).  The  necessary  conditions  for  this  coarse  grid  problem  are 

AcXc  -  Bcuc  =  g% 

A\f  +  CtCcx*  =  g%  (4.10) 

B'PC  =  gZ 


(4.8) 

(4.9) 
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which  are  the  FAS  equations  for  the  fine  grid  necessary  conditions.  In  particular  they 
have  the  property  of  introducing  no  change  to  the  fine  grid  problem  once  the  minimum 
has  been  reached. 

Once  an  approximation  to  the  coarse  grid  problem  is  found  it  is  used  to  correct  the 
fine  grid  approximation  by 


x*  <r-  x*  +  Px(xc  -  Rxxf) 

ui  <_  uf  4.  Pu(uc  -  RuU*)  (4.11) 

pf  <_  pi  +  Pp(pc  _  Rppf) 

4.2  Relaxation 

The  general  equations  for  state  and  costate  variables  on  any  of  the  grids  involved  in 
the  process  have  the  form, 


Ax  -  Bu  =  fx 
A*p  +  C*Cx  =  /p 

These  two  equations  are  being  relaxed  on  all  levels  using  one  of  the  standard  relax¬ 
ation  methods  for  elliptic  problems.  That  is,  Gauss-Seidel  or  damped  Jacobi  in  case  A 
is  a  discretization  of  a  scalar  PDE,  and  Collective  Gauss-Seidel  (CGS)  or  Distributed 
Gauss-Seidel  (DGS)  when  systems  of  elliptic  PDE’s  are  involved.  So  the  problem  of 
relaxation  for  x  and  p  is  quite  standard.  An  important  property  of  such  relaxation 
schemes  is  that  for  a  wide  class  of  discretization  methods  (  h-elliptic  discretizations) 
they  have  fast  convergence  for  h-non-smooth  errors  (defined  later)  (see  [B2]).  This 
implies  that  the  errors  in  the  x  and  p  equations  will  have  smooth  residuals  after  just  a 
few  relaxation  sweeps.  This  does  not  mean  that  the  errors  for  the  optimization  prob¬ 
lem  are  smooth.  The  relaxation  of  equations  (4.12)  serves  as  a  process  that  brings  the 
solution  closer  to  the  constraint  surface  defined  by  the  first  equation  there. 

As  far  as  updating  the  control  variables  the  situation  is  more  complicated.  Note 
that  a  small  change  in  any  of  the  control  unknowns  Uj  implies  a  change  in  x  and  p 
everywhere  in  the  domain.  Good  efficiency  for  the  full  algorithm  is  to  be  obtained  if  the 
changes  in  *  and  p  that  result  from  changes  in  u,  can  be  calculated  with  a  computational 
cost  which  is  not  much  greater  than  that  of  relaxing  the  equations  (4.12)  once. 

The  smoothness  assumption  on  the  coefficients  of  the  operator  A  implies  that  non¬ 
smooth  errors  correspond  to  non-smooth  residuals  and  vice  versa.  Non-smooth  errors 
are  also  fast  to  converge,  for  h-elliptic  discretizations,  as  mentioned  before  (this  is  not 
the  case  for  quasi-elliptic  discretization  schemes).  This  results  in  the  following  basic 
rule  for  updating  the  control  variables. 

•  Make  changes  in  the  control  variables  that  result  in  non-smooth  residuals 

Note  that  a  change  in  any  of  the  control  variables  may  result  in  smooth  changes 
on  one  grid  but  non-smooth  changes  on  a  coarser  grid.  This  suggest  that  the  different 
control  variables  will  be  relaxed  on  possibly  different  grids  depending  on  their  scale  of 
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influence  in  the  solution.  To  quantify  more  precisely  these  ideas  we  need  the  following 
definitions. 

Definition:  A  function  if>h(0  will  be  called  h-non~smooth  at  {  €  Slh  iff 

I (f*hTj>h)(0  ~  1  >  a#fc(OI  0  <  a0  <  a  <  l,a0  =  0(1),  (4.13) 


where  is  a  local  averaging  operator. 

Definition:  A  function  ij>h(£)  will  be  called  uniformly  h-non-smooth  in  SI  iff  it  is  h- 
non-smooth  in  every  point  £  £  Slh 

Thus  if  i/itf  is  uniformly  h-non-smooth  its  effect  on  the  solution  x  and  p  can  be 
calculated  accurately  with  just  a  few  relaxation  sweeps.  Typically  2-3  relaxation  sweeps 
are  enough  to  reduce  the  errors  by  an  order  of  magnitude.  This  will  be  the  key  to  our 
method  of  updating  the  control  variables  Uj.  Let 

S*  =  €  Slh  :  il>j  is  h-non-smooth  at  (4-14) 

If  h  correspond  to  the  coarsest  level  and  then  we  define  Ej  =  Slh.  Define 


dkk(0 


Dfc(0  =  diag(41(0,...,4^)) 

if  ip1?  is  h-non-smooth  at  £ 


={: 


0  otherwise 


(4.15) 


The  qxq  matrix  Dh  will  be  used  in  defining  the  relaxation  for  the  control  variables. 

In  case  the  functions  ifij  are  uniformly  h-non-smooth  the  update  of  the  control 
variables  becomes  much  simpler.  The  perturbation  in  u  that  is  done  in  minimizing 
the  functional  is  by  jDh(Bh)*ph,  followed  by  a  few  relaxation  sweeps  of  the  state  and 
costate  equation  in  a  vicinity  of  (J  yielding  a  good  enough  approximation  to  xh 
and  ft1.  The  actual  perturbation  is  then  chosen  so  as  to  minimize  the  functional  or  the 
other  choice  as  explained  in  section  2.  Note  that  the  multiplication  by  D*1  causes  the 
perturbation  in  the  right  hand  side  to  be  h-non-smooth  on  the  given  level.  That  is, 
each  control  variable  is  being  updated  on  the  appropriate  grid. 

For  il>j  with  a  singularity  a  slightly  different  approach  is  taken.  Here  the  effect  on  the 
solution  is  smooth  away  from  the  vicinity  of  non-smoothness  of  rj>j.  This  follows  form 
basic  theories  of  elliptic  partial  differential  equations  [ADN].  Thus,  if  the  observation 
operator  C  is  supported  far  enough  from  the  singularities  of  V’j,  even  coarse  grids  can 
compute  accurately  the  effect  on  Cx,  which  is  the  important  quantity  in  updating 
the  Uj.  However,  when  the  singularities  are  close  to  the  support  of  C  this  is  not  the 
case  and  some  refinement  has  to  be  done  locally  in  order  to  account  for  the  local  non¬ 
smoothness.  In  practice,  a  few  points  around  the  singularity  need  to  be  relaxed,  further 
points  are  relaxed  on  coarser  grids.  The  finest  grid  which  needs  to  be  involved  in  the 
process  with  a  given  singularity  is  determined  from  the  accuracy  achieved  in  Cx  which 
can  be  estimated  using  the  quantities 

|| Ch$h  -  Chi2h ||  (4.16) 


and  the  knowledge  about  the  accuracy  of  the  discretization  involved.  The  later,  if  not 
known,  can  be  computed  from 

lie'1*'1  -  (4  17n 

\\Chx2h  -  ChiAh\\  K  ' 
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In  any  case,  the  work  performed  in  the  calculation  of  xk,ffr  should  not  involve  more 
than  one  multigrid  cycle,  done  as  an  FMG  cycle  used  with  local  refinements.  The  exact 
formulation  of  the  algorithm  for  updating  the  control  variables  in  the  most  general  case 
will  be  presented  elsewhere.  We  denote  the  process  of  updating  the  control  variables 
in  the  direction  u  by 


(x,p,u)  <—  MinRel(x,p,u,u)  (4.18) 

Having  defined  the  relaxation  and  the  coarsening  process  the  full  algorithm  can 
now  be  described. 


4.3  Unistep  Algorithms 

Let  Vk,Wk,Uk  and  Hk  be  Hilbert  spaces  approximating  V,  W,W  and  H,  respectively. 
Let  Ak,Bk  and  Ck  be  operators  as  follows  Ak  :  V*  W*,  Bk  :Uk  -*  W*,  Ck  :  V*  -» 
Hk,  approximating  the  operators  A,B,C,  respectively.  Assume  also  that  interpolation 
and  restriction  operators  Pk,Pk,Pk,Rk,Rk,R!kl,Rk,Rk,Rk  are  given  as  Rk  :  V*+i  -*■ 
V*,  12*  :  Wfc+1  -  Wfc,  Rku  :  Uk+i  -  Uk,  Rk  :  Vfc+1  -  V*.  Rk  :  Wfc+i  -»  Wfc, 
Rt  •'  *4+1  -*  *4,  PX  ■  Vk  -»  V*+1,  Pp  :  Wk  -»  W/fe+1,  :  *4  -»•  *4+1-  Also  assume  the 
the  qxq  matrices  Dk  (defined  in  section  4.2)  axe  given. 

On  all  levels  minimization  problems  as  follows  are  given, 


<  Ckxk,xk  >  -2  <  ®*,s*  >  -2  <  uk,gk  > 
Akxk  =  Bkuk  +  gk 


(4.19) 


Observe  that  this  is  equivalent  to  having  the  following  systems  of  equations  on  all 
levels 


Akxk  =  Bkuk  +  gk 

A*kpk  +  C*kCkxk  =  gk  (4.20) 

Btf  =  gk. 

Let  (®*,p*,rt*)  be  an  approximate  solution  of  (4.19).  We  define  next  an  algorithm 
for  improving  it,  denoted  by 


(«*,?*,«*) «-  MG{xk^tuk,gkx,gk,gk)  (4.21) 

If  k  =  1  Then 

1.  relax  the  first  two  equation  in  (4.20)  until  convergence. 

2.  iterate  until  convergence: 

2a.  perform  the  cycle 

(®*»p\«*)  <-  Afiniiei(®fc,pfc,ttfc,2?*(p*-HJpfe)) 
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2b.  relax  the  first  two  equations  in  (4.20). 

Else 

1.  perform  the  following  v\  times 

la.  relax  the  first  two  equation  in  (4.20) 

lb.  perform  the  cycle 

(**,?*,«*)  MinRel(xkipft,uk,Dk{gk-Blpk)) 

2 .  Let  k  =  k  -  1 ,  and 

&  =  Rkx  (gk+1  -  Ak+1xk+' +  Bk+m*1)  +  AkRkxk+1 
9kP  =  Rkp  (9k+1  ~  A*fc+1p*+1  +  Ci+1Ck+1uk+l)  +  AkRk p*+* 

9t  =  K  (9t+1  ~  +BkRkPk+1 

3.  perform  7  times  the  cycle 

(xk,pk,uk)  <-  MG{xk,j^,uk,fkJk,fk) 

4.  correct  fine  grid  solutions 

j.fc+1  <_  xk+ 1  +  pk(xk  _  Rk  j.fc+1) 

p*+ 1  <-  pk+1  +  Pk(pk  -  Rk  pfc+1) 
uk+ 1  4-  uk+1  +  Pk(uk  -  Rk  u*+1) 

5.  perform  the  following  v\  times 

5a.  relax  the  first  two  equation  in  (4.20) 

5b.  perform  the  cycle 

(xk,pk,v.k)  <-  MinRel{xk,jft,uk>Dk{gk-Blpk)) 

End 

In  order  to  obtain  full  efficiency  the  algorithm  starts  at  the  coarsest  level,  where 
each  time  a  refinement  is  done  followed  by  a  fixed  number  of  MG  cycles.  This  is  the 
N-FMG  algorithm  which  is  defined  next. 

1.  Solve  (4.19)  for  *  =  1 

2.  *  =  *  +  1,  pfc  =  n*-V-1,u*  =  p*-1u*-1. 

3.  Define 

4.  Perform  N  time  the  cycle 
(xk,p^,uk)  *-  MG(xk,pfe,uk,fk,fktfk) 

5.  If  k  =  M  stop,  else  goto  2. 

At  the  end  of  this  algorithm  an  approximate  solution  to  the  minimization  problem 
on  level  k  —  M  is  given. 

5  Numerical  Examples 

Numerical  experiments  were  conducted  with  scalar  elliptic  problems  governed  by  the 
Lapladan  in  the  unit  square  with  Dirichlet  boundary  conditions.  That  is,  using  the 
notation  of  section  2, 


12 


(5.1) 


fi  =  {(*>2/) ;  o  <  v>y  <  1} 

V  =  w  =  {<f>e  L2(Q) :  tf|an  =  0} 

A  =  A 

H  =  L2(T)  r  =  {(x,0):0<s<  1} 

C  =  d/dn\T 

The  control  space  was  Mq.  Different  cases  were  considered  for  (ipj  j  =  1,. . .  ,q)  , 
in  the  definition  of  B. 

Uniform  grids  were  used  in  the  discretization  with  the  standard  5-point  formula  for 
the  Laplacian,  namely, 


Slh  =  {( i/NJ/N )  :  0  <  ij  <  N} 

Th  =  {(i/N,0):0<i<N}  (5.2) 

(Akxh) u  =  MxUi  +  x!*i.i  +  Xh-1  +  xti+i  -  4Xb) 

where  Nh  =  1,  h  being  the  discretization  parameter.  The  observation  operator  was 
discretized  using  a  3-point  formula, 


(chx% o  =  i(l-5Xf0  -  2X*  +  .5X1:,) 


(5.3) 


and  the  functions  were  discretized  by  simple  injection. 

With  this  discretization  the  discrete  solution  is  expected  to  have  second  order  ae 
curacy.  That  is, 


\\Xh-X\\<fixh2 
||«*  -  «||  <  Puh2 


(5.4) 


where  PX,PU  are  constant  independent  of  h  (they  depend  on  high  order  derivatives  of 
the  state  variable).  It  is  enough  to  solve  the  discrete  problem  to  an  accuracy  satisfying 


II**  -  **||  «  ||**  -  *|| 

||*fc  -  uh\\  <  ||u*  -  tt|| 


(5.5) 


Here  the  tilde  quantities  represent  the  current  numerical  approximations  to  the  solution 
of  the  discretized  problem.  The  quantities  |j**  -  *||,  ||u*  -  u||  are  called  discretization 
errors. 

In  all  the  examples  reported  below  the  algorithm  described  in  section  5  was  used 
with  the  following  parameters:  v\  =  2, 1/3  =  1.  The  restriction  operators  for  both  the 
residuals  and  the  full  solutions  Rx,  R*,  **,  R*  were  the  9-point  full  weighting,  that  is 
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1  2  1 
2  4  2 
1  2  1 


(5.6) 
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The  interpolation  operators  P*+1 ,  P*+1  were  the  standard  bi-linear  interpolation,  II* ,  II* 
the  bi-cubic  interpolation  and  -R*  =  Ru  =  I ,  the  identity  operator.  The  mesh  size  on 
level  k  was  2"(*+1\  k=l  being  the  coarsest  level.  Results  for  3-FMG-V(2,l)  cycle  are 
given.  In  each  table  results  for  the  L3  residuals  rx,  rp,  ru  of  the  state,  costate  and  the 
control  equations  are  given.  Their  exact  definition  is  given  by 

Ml  =  £  ((jl)«  +  (^«*)«  -  (5.7) 

«J=1 

Ml  =  £  ((st)ij  -  (<7{C»x* (5.8) 

Ml  =  £((*J).  -  (VP*),)2  (5.9) 

1=1 

where  the  scaling  of  h\  in  these  definitions  is  used  so  that  these  norms  approximate 
the  continuous  norms  and  so  residuals  on  different  levels  can  be  compared. 


1.  Distributed  Control 

In  the  distributed  control  case  we  have  worked  with  the  following  problem 

“S*  fr(jh  -  d)2da  (51°) 

Ax  =  £<=1  Uiipi 

x|an  =  0.  (5.11) 

ipj  =  ain(njir^i)sin(mjtr^2)  nj>mj  integers, 

Table  1  shows  the  results  for  a  one-dimensional  control  case,  that  is,  q  =  1  where 
mi  =  Tix  =  1,  i.e.  a  smooth  case.  This  is  the  simplest  of  all  cases  and  is  given  here 
mainly  for  reference.  Residuals  of  the  state,  costate  and  control  equations  are  given 
separately.  The  error  in  the  control  (relative  to  the  true  differential  error)  is  given  as 
well.  Observe  that  a  1-FMG-V  cycle  gives  solutions  to  the  levels  of  discretization  errors 
on  fine  levels,  as  ||v  -  tie*oct||  reaches  its  minimum  in  essentially  one  cycle.  The  0(h2) 
behavior  of  the  error  is  clear  from  the  results,  reflecting  the  order  of  the  scheme  used. 

Table  2  shows  the  results  of  a  similar  experiment  in  which  the  control  is  three 
dimensional  (q  =  3)  with  n\  =  ?»3  =  l,mi  =  mj  =  1,  rtj  =  m3  =  2.  The  behavior  here 
is  similar  to  the  previous  case.  The  0(h2)  is  evident  also  here. 


2.  Boundary  Control 

The  boundary  control  case  problem  was 


(5.12) 
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(5.13) 


Az  =  0 
*|en  =  E?=i  «• 

i>i  =  a*n(n,ir£i)|ao  integer, 

Tables  3  and  4  shows  results  for  this  problem.  The  first  with  q  =  l,ni  =  1,  the 
other  with  q  =  2,ni  =  l,nj  =  2.  Also  in  this  case  the  results  are  basically  optimal,  a 
solution  to  the  levels  of  discretization  errors  are  obtained  in  1-FMG-V(2,1).  The  0(h2) 
convergence  toward  the  differential  solution  is  evident  in  both  cases.  Observe  that  the 
initial  residuals  on  finer  levels  start  are  smaller  by  a  factor  which  is  close  to  four.  This 
is  typical  to  examples  with  smooth  solution  and  proper  FMG  interpolation  operators 
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S.  Pointwise  Control 

The  pointwise  control  case  problem  was 


(5.14) 


Az  =  £J=1  unl>i 

z|ao  =  0.  (5.15) 

=  *(*-&).&€  n 

Tables  5,6  and  7  show  results  for  non-smooth  control.  Results  are  given  for  one 
two  and  three  delta  functions,  i.e.,  q  =  1,2,3.  The  location  of  the  delta  function  is 
given  in  each  table.  Here  the  solution  is  less  smooth  than  before.  Still  the  results 
show  essentially  the  same  behavior  as  in  the  smooth  case.  Observe  that  the  initial 
residuals  on  the  different  levels  are  of  the  same  order,  reflecting  the  non-smoothness  of 
the  solution.  This  behavior  cannot  be  improved  by  using  a  higher  interpolation  in  the 
refinement  stage  of  the  FMG  algorithm,  although  local  relaxation  in  the  vicinity  of  the 
singularities  can  improve  the  results  we  have  not  experimented  with  such  ideas  here. 

In  some  of  the  experiments  it  can  be  observed  that  one  of  the  three  residuals  shown 
is  increased  in  some  of  the  cycles.  This  is  reasonable  since  the  three  residuals  are 
coupled  and  it  is  only  their  sum  which  goes  down  in  each  cycle.  The  results  presented 
by  all  tables  clearly  demonstrate  the  effectiveness  of  the  method  developed  here  which 
aimed  at  the  full  optimization  problem,  therefore  leading  to  one-shot  solution  of  these 
problems.  1-FMG-V  cycle  is  basically  enough  to  reach  below  the  levels  of  discretization 
errors. 
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level  |  cycle  no. 


5 


2 


|u  —  uexact 


.8166 


.589e-01 
.612e-01 
.568e— 01 
.581e-01 
.576e-01 


.174e— 01 
.149e— 01 
.147e-01 


.329e— 02 
.305e— 02 
.301e-02 


.337e-06 


.591e— 02 
.104e-03 
.216e-03 
.592e-04 
.250e-04 


.860e-03 

.648e-03 

.121e-04 


.125e— 03 
.324e-04 
.407e— 04 


.121e-05 


.121e+00 

.609c-02 

.373e-02 

.996e-03 

.446e-03 


,522e-01 

.292c-02 

.136e-03 


.357c-01 
.132e— 02 
.683c-04 


.397e-09 


.415e— 03 
.975e-04 
,485c- 04 
.183c-04 
.737e-05 


.319c-04 

.209c-06 

.7556-07 


.434e-05 
.371e— 06 
.294e— 07 


Table  1:  q  =  !,«!  =  l,mj  =  1 


level  cycle  no. 


5 


2 


ti  —  ttexoct 


.246e— 06 


.8876-02 

.3786-02 

.169e-02 

.2376-03 

.905e-04 


.246e-02 
.418e— 03 
.113e-03 


•380e-03 

.595c-04 

.383e-04 


.110e-05 


.207c+00 
.592e-01 
.258c-01 
.359e— 02 
.130e-02 


.110e+00 

114e-01 

.294e-02 


,618e-01 

.316e-02 

.167e-03 


.374e-08 


.654c-03 

.200c-03 

.103e-03 

.271e-04 

.908e-05 


4 

4 

5 


.835e-05 

.110e-05 

.166e-06 


.2536+00 


.174e+00 

.117e+00 

.994e-01 

.100e+00 

.102e+00 


.427e-01 

.317e-01 

.293c-01 


.8586-02 

.630e-02 

.618e-02 


Table  2:  g  =  3,nj  =713  =  1,  mi  =  m3  =  1,  nj  =  m3  =  2 


Table  3:  q  =  l,m  =  1 


level  cyde  no. 


5 


2 


u  —  Vgxact 


JB 


.631e-06 

.800e-06 

.296c-13 

•334e+00 

.354e+00 

•786e+01 

.288e-01 

.149e— 01 

.206e+00 

.388e+01 

.896e-02 

.150e+00 

.120e+00 

.229e+01 

.307e-02 

.754e— 01 

•706e-01 

.134e+01 

.106e-02 

.117e+00 

.415e-01 

•792e+00 

.367e-03 

.9306-01 

.523e-01 

•250e+01 

.9246-03 

.3926-01 

.106e-01 

.526c+00 

.4466-04 

.2986-  01 

.212e-01 

.104e+00 

.1686-05 

.2816-01 

.120e-01 

•126e+01 

.394e-05 

.5286-02 

.628e-03 

.668e-01 

.9586-08 

.600e-02 

.274e-03 

.466e— 02 

.313e-10 

.594e-02 

Table  4:  q  =  2,r»i  =  l,n2  =  2. 


level 

cyde  no. 

mmm 

■n 

ru  3 

U  —  U exact  2 

1 

5 

.211e— 05 

.921e-05 

.347e-06 

.8646+00 

2 

.619e+00 

.437e+01 

.363e-01 

.390e+00 

.2366+00 

.223e+01 

.256e-01 

.1816+00 

.147e+00 

.144e+01 

.166e-01 

.5366-01 

.927e-01 

.921e+00 

.105e-01 

.4256-01 

5 

.598e— 01 

•591e+00 

.677e-02 

.9096-01 

3 

1 

.391e+00 

.389e+01 

.1486-01 

.537e-01 

2 

.389e-01 

.200e+00 

.160e-02 

.3736-01 

3 

.138e-01 

.336e-01 

.523e-03 

.3126-01 

mm 

1 

.660e+00 

.3616+00 

.5876-03 

.6736-02 

1  1 

2 

.168e-01 

•187e-01 

.405e-04 

.5116-02 

1 

3 

.708e-03 

.120e-02 

.116e-06 

.5096-02 

Table  6:  q  =  2, =  (.75,  .5),  6  =  (.5,  .5) 
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level  cycle  no. 


IMa 


.131e— 05 


.816e+00 

.293e+00 

.215e+00 

.131e+00 

.110e+00 


.396e+00 
.431e-01 
.244e— 01 


,683e+00 

.402e-01 

.452e-02 


TA  —  Inexact 


19 


.953e-05 


.447e+01 

.245e+01 

.177e+01 

.119c+01 

.883c+00 


.412e+01 

.214e+00 

.4766-01 


.3736+00 

.341e-01 

.533e-02 


.334e— 05 


.477e— 01 
.319e-01 
.243e— 02 
.162e-01 
.131c— 01 


.155e-01 
.236e— 02 
.103e-02 


.132e-02 
.258c-03 
.309e— 04 


.115e+01 


.531e+00 
.273e+00 
.848e— 01 
.345e— 01 
.127e+00 


.793e— 01 
.640e-01 
.531e-01 


.168c-01 
.880c-02 
.788e— 02 


Table  7:  q  =  3,6  =  (.75,  .5),  6  =  (-5,  .5)6  =  (-25,  .5) 
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